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We present a self-contained description of everything needed to write a program that calculates 
the CMB power spectrum for the standard model of cosmology (ACDM). This includes the equa- 
tions used, assumptions and approximations imposed on their solutions, and most importantly the 
algorithms and programming tricks needed to make the code actually work. The resulting program 
is compared to CMBFAST and typically agrees to within 0.1% - 0.4%. It includes both helium, 
reionization, neutrinos and the polarization power spectrum. The methods presented here could 
serve as a starting point for people wanting to write their own CMB program from scratch, for in- 
stance to look at more exotic cosmological models where CMBFAST or the other standard programs 
can't be used directly. 



I. INTRODUCTION 

Since the first detection of anisotropies in the Cos- 
mic Microwave Background (CMB) by the Cosmic Back- 
ground Explorer (COBE) satelhte in the early 90's 0, 
there has been considerable activity in this field of cos- 
mology throughout the world. With the more accurate 
measurements of the Wilkinson Microwave Anisotrow 
Probe (WMAP) 0,13 and the future Planck satellite [J, 
we are entering the era of precision cosmology. Results 
from different kinds of experiments seem to converge to 
what is being referred to as the Standard Model of Cos- 
mology 0, describing the history of the universe from 
inflation through Big Bang Nucleosynthesis to the release 
of the CMB radiation during recombination. With the 
growing precision of the observational data comes also 
the need for fast and accurate theoretical calculations of 
CMB power spectra. Often one seeks the best fit to ob- 
servations of a model with several parameters, requiring 
typically hundreds of spectra to be calculated and com- 
pared to the data. 

Of course, several standard computer programs that 
calculate CMB power spectra are already available. 
The most commonly used include CMBFAST H 0, 
CMBEASY 0, and CAMB 0. These are all excellent 
programs when calculating power spectra for the ACDM 
model, including also some extensions like a quintessence 
field, hot dark matter or a simple fifth dimension. One 
may therefore wonder what the point is writing a new 
program from scratch, a task which obviously requires 
quite a lot of work. The need may arise when considering 
more exotic cosmological models which are not included 
in any of the standard programs. This could for instance 
include quintessence with a non-trivial coupling to other 
fields, extra dimensions with non-trivial geometry, or a 
model with varying constants of nature. In these cases 
one has the choice of either extending existing code, or 
writing a completely new program. The former can cer- 
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tainly be the easiest solution in some situations, but not 
necessarily all. The problem with updating existing code 
is that one probably doesn't know exactly how it works, 
and it is therefore difficult to make changes without doing 
something wrong. By instead writing new code, starting 
with a simple model and comparing the results to exist- 
ing programs, it will be much easier to later extend the 
program since one then knows precisely what every line 
of code is doing. By knowing the code in detail, more 
confidence can also be put in the results. Obviously a lot 
of work is needed to write a CMB program from scratch. 
Fortunately this work is not a complete waste of time, 
since a lot of new insight about the underlying physics 
can be obtained in the process. 

The purpose of this text is to provide a self-contained 
collection of all the ingredients needed to write a pro- 
gram that calculates the CMB power spectrum. This 
includes both the equations governing the physics, any 
assumptions or approximations used in their solutions, 
and pointing out what algorithms and tricks to use when 
implementing the equations in a computer program. The 
main focus will be the practical computer implementa- 
tion of the equations, not their derivations, for which we 
instead point the reader to the references. Only new or 
less readily available derivations will be included. We are 
only assuming that the reader has some basic knowledge 
of CMB physics, and some experience with a high-level 
programming language ^. Most of the presentation will 
follow the notation and conventions of Dodelson . 

In section^lwe start by going through the background 
cosmology, which includes background geometry and the 
recombination history of the universe, and in section Hill 
we introduce perturbations to this background. The 
equations of motion for the perturbations are given by 
the various Boltzmann equations. We then state the ini- 



^ We will not show any code written in a spesific programming 
language, only the general alg orit hms used. For the record, we 
have used Delphi for Windows llj in this work, but C or Fortran 
(or similar languages) are equally well suited. 
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tial conditions used, and the approximations used dur- 
ing tight coupling. In section IIVI we go from the per- 
turbations of the CMB temperature to the spectrum of 
C/'s. The most important programming techniques are 
mentioned in section including the cutoff scheme for 
the Boltzmann hierarchies, how to integrate the various 
equations numerically, and the normalization of the spec- 
trum. The resulting spectra are compared to CMBFAST 
in section IVII In section IVIII the program is extended 
to include a few more effects, including helium, a simple 
model of reionization, massless neutrinos, and the (E- 
mode) polarization power spectrum. Finally we conclude 
m section rvml 



II. BACKGROUND COSMOLOGY 



The various "time variables" t, rj, a and x are related 
through the useful equations 



dt dx ^ 

dj] ' dt ^ drj ' 

d ^ d d ^ d 
dt dx ' dr/ dx 



(4) 



We will also need an expression for the conformal time 
as a function of the scale factor in our calculations: 



via) 



da' 



a'nia') 
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This integral is easily calculated numerically. Note that 
a7i{a) HqV^ as a 0, so there's no problem with 
convergence. 



Before looking at perturbations, we must determine 
the background. This consists of two parts: The easi- 
est part is the background geometry, which is given by 
the standard Friedmann-Robertson- Walker (FRW) met- 
ric. The more difficult part is the recombination history 
of the universe, which involves finding the number of free 
electrons and electrons bound to neutral atoms as a func- 
tion of time. We need this to determine the coupling 
between photons and baryons. 



A. Background geometry 



The background geometry is given by the FRW metric 



ds^ = -dt^ + a^{t)6ijdx'dx^ 

= a^{r]) {-drf' + Sijdx^dx^) , 



(1) 



where t is the physical time and rj the conformal time, 
and a is the scale factor describing the expansion of the 
universe, which we assume is spatially flat (fc = 0). The 
expansion is given by Friedmann's equation 



1 da 



H = -— = Hoy^{nrn + flb)a-^ + flra-^ + nA , 
a dt 

n^^-^^^=aH 
a dr] a 



= Ho^y{n^ + r2b)a-i + f2r-a-2 + f^Afl^ , (2) 

where the dot means the derivative with respect to con- 
formal time, and we assume that the universe consists of 
cold dark matter (CDM, m), baryons (6), radiation (r), 
and a cosmological constant (A) . Hq is the current value 
of the Hubble constant. We also introduce the logarithm 
of the scale factor, 



In a , 



d 

dx 



(3) 



B. Recombination 

In the early universe all atoms were fully ionized, giv- 
ing a strong coupling between the baryon and photon 
plasma due to Thomson scattering. When the tempera- 
ture dropped below ^ 3000 K neutral atoms were formed, 
and the universe became transparent. The CMB pho- 
tons we observe today have travelled more or less freely 
through the universe since they were last scattered during 
recombination. The optical depth r back to conformal 
time rj is given by 

rvo 

t{vi) = / necrradr]' , 

T = -Heara ^ r = — — , (6) 

n 

where rif, is the number density of free electrons, 



6.652462 x IQ-^s 



(7) 



the Thomson cross section, and 770 the conformal time 
today, r/o — rj{a = 1). We define the visibility function 



g[x)^-r'e-^^^^. (8) 



n{x) 

The visibility function is normalized as 



g{vi)dr] = / g{x)dx = 1 



(9) 



and can therefore be interpreted as a probability distri- 
bution, namely the probability that a CMB photon ob- 
served today was last scattered at conformal time 77. The 
function g has a relatively sharp peak at a certain red- 
shift, of order z ^ 1100, which we therefore call the time 
of recombination. Most CMB photons were last scattered 
around this time. 
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The difficult task is to calculate the electron density 
Up. We define the free electron fraction 



n,. 



rib 



(10) 



where the total number density of hydrogen, riu , is equal 
to the baryon number density Ub when we ignore helium 
(see section IVlI A\i . Ignoring also the small mass differ- 
ence between free protons and neutral hydrogen, we have 



riH = rib 



rtbPc 



Pb 



Pc 



3Hl 
8ttG 



(11) 



Here mpj is the mass of the hydrogen atom, and pc the 
critical density today. At early times all hydrogen is com- 



pletely ionized, so ~ 1 and rig 



whereas at late 



times Xe ^ 1 (but does not approach zero). 

Before recombination the electron fraction can be ap- 
proximated by the Saha equation 0, [13. IT^ : 



X! 



l-Xp 



J_ f meTb 

Ub 



27r 



3/2 



(12) 



Here Tb is the baryon temperature, and eo — 
13.605698 eV the ionization energy of hydrogen. Dur- 
ing and after recombination, however, we must use the 
more accurate Peebles equation [ill uM 



dXe 
dx 

where 



CrjTb) 

H 



(3{Tb){l-Xe)-nua^^\Tb)X'^ , (13) 



a(Tfc) 



A, 



A, 



A 



Ao^+p(^)iTb) ' 



.227 s~^ , 



A. 



H- 



(3eo)^ 
(87r)2n 



Is 



~ (1 - Xe)nn , (Tb) = P{Tb)^'^>'^^^ , 



P{Tb) = a'^^\Tb) 



a^'\Tb) 




02 (Tfc) - 0.448 ln(eo/r6 



(14) 



Naively, one would probably think that recombination 
occurs when T ~ eg ~ 157 900 K when looking at p2|l or 
(tO|l . However, it is delayed until T ~ 3000 K because of 
the large photon to baryon number ratio. 

It is difficult to integrate the Peebles equation l|13() 
numerically at very early times, but this is also the place 
where the Saha equation (|12(l is a good approximation. In 
the numerical calculation we therefore use Saha until the 
electron fraction X^ has been reduced to, say, 0.99, and 
then switch to Peebles using Saha as the initial condition. 
Figure Q shows the numerical solution X^, as a function 
of redshift z — a"^ — 1, and figure|21the resulting optical 
depth and visibility function, all for the model 

Hubble constant: h = 0.7 , 

iHo = h- 100 kms-i Mpc-i) 
CMB temperature: Tq = 2.725 K , 

{=> a- = 5.042 X IQ-^) 
Baryon density: ilf, 
CDM density: n 
Spectral index: 
Helium mass fraction: Yp 
Vacuum density: J^a 



0.046 , 
= 0.224, 

= 1, 
0, 

= i-{nr- 

= 0.72995 , 



(16) 



which we call the "default model" for the rest of this text. 
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FIG. 1: The free electron fraction as a function of redshift, 
using the Saha approximation 11211 until z — 1587.4 where 
Xe = 0.99, and then integrating the Peebles equation 11311 . 



The various terms here are described in more detail in 
|12|. The baryon temperature Tb has a non-trivial time 
evolution, and is g iven by a differential equation which 
couples to Xe [131 • Thus, we actually have a compli- 
cated coupled system of differential equations for both 
Xe and Tb- However, the error when setting the baryon 
temperature equal to the photon temperature through- 
out recombination turns out to be only of order 10~^ 
We therefore use the approximation 



To 
a 



Tq = 2.725 K . 



(15) 



III. PERTURBATIONS 

A. Definitions 

Having determined the background cosmology, we can 
now turn to the perturbations. We use the Newtonian 
gauge and write the perturbed metric as 



_ / -(1 + 2^-) 
^''""l, a2<5y(l + 2$) 



(17) 
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FIG. 2: The top figure shows the optical depth r (solid line) 
and |r'j (dashed line) as functions of x. The bottom fig- 
ure shows the visibility function g — — r'e"'^ (solid line), 
its derivative g'/W (dashed line), and its double derivative 
g" /300 (dotted line), also as functions of x. The scaling is 
chosen to make the curves fit into the same figure. The vis- 
ibility function has a peak at a; = —6.984, corresponding to 
a redshift z — 1078. The peak of g has both a finite width 
and a clear asymmetry, and is therefore not particularly well 
approximated by a delta function. 



B. Perturbation equations 

The equations of motion for the perturbations fol- 
low from the Boltzmann equations for photons, CDM, 
baryons and neutrinos. Since CDM and baryons are non- 
relativistic, we can take the first two moments of their 
equations instead of keeping an arbitrary direction de- 
pendence, obtaining equations for the density and veloc- 
ity. In addition to the Boltzmann equations, Einsteins 
equation gives two equations for the two gravitational 
potentials. In total, the evolution of the perturbations is 
governed by the following system of equations : 



9 -I- ikfiQ = -$ - ikfi'i' 

— f 9o — 6 + ifJ-Vb 
1 



Qp + ik^Qp = — f 



:{1-V2)U 



6 — kv 
Sb - kvb 
Vb + Hvb 
jif + ik^Af 



Hv 



-3$, 
-3$, 

-k^ + rR {vb + iQi) 
-<i) - ikii'^ , 



47rG'a^ 



pS + pbSb 

+ 4p^eo4 



n = 62 



e 



-327rGa2 
R 



4pr 

3/Ob 



PrQ2 + PuJ^2 



iVLba 



(20) 



We are therefore only considering scalar perturbations. 
Perturbations to the photons are defined as the relative 
variation of the photon temperature: 



T(fc,/i,r;) =T(")(77) [l + e(fc,M,7?) 



■P 



kp 



(18) 



We will be working in Fourier space throughout this text, 
with k as the Fourier transformed variable of the position 
X. The momentum of the photon itself is p. Note that the 
perturbation Q depends only on the direction of p, not 
its magnitude. The direction dependence is what leads to 
anisotropies in the CMB field. The photon perturbation 
is expanded in multipoles: 



°° 2/ + 1 



(19) 



1=0 



where Vi{p) are the Legendre polynomials. In addition to 
the temperature perturbation, there's also perturbations 
to the photon polarization, which we denote by Qp{p). 
(See section FVlI Dl for more details on polarization.) 



Here S and v are the density perturbation and velocity of 
CDM, 6b and Vb the same for baryons, and Af (massless) 
neutrino perturbations. Compared to 0, eqs. 4.100 - 
4.107] we have defined R — > 1/R, and v —^ iv, Vb — > ivb 
to make the velocities real. Expanding in multipoles, the 
equations for O, Qp and J\f turn into the hierarchies 



eo + fcei = 



Ik 



21 + 1 



&1-1 



{l + l)k 
21 + 1 



e 



i+i 



Ik 



21 + 1 



e 



i-i 



il + l)k 
21 + 1 



Ml - ^Mo + yAAa 



Ik 



21 + 1 



/-I 



(/ + l)fc 



21 + 1 



3 



il>2) 



Ql + ^Vb 



10 



-US, 



1,2 



-n 



er -—ns, 



10 



1,2 



a>i) 

N'o + kNi = -$ , 
k 
3 





(/ > 2) (21) 
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Using X as the time variable and rearranging the equa- 
tions, we obtain their final form: 



0; 
Q'pl 



n 



61 - , 



3n ' m 



Ik 



i2i + i)n 
k 



p 



ei-i 
1 



k T 

* + T 

{l + l)k 

{2i + i)n 



3 " 



e 



10 



1,2 



l>2, 



n 



Ik 



{2i + i)n 



{l + l)k 

{2i + i)n 



er - — n^, 



10 



1,2 



M2 + — ^' 



' {2i + i)n 



m 



m 
(/ + i)fc 

(2Z + 1)W 



l>2. 



v' = -v--^, 



H 



-Vb 



^ T 
^ 

n 



r'i? (391 + i;^) , 



Hi ^ 



2W 



k'^a?' 



(22) 



The expression for 5" is just an algebraic equation, so 
this expression should simply be inserted into all the 
other equations when needed. Also, the expression for 
should be calculated first and used in all the other equa- 
tions, so that we obtain a system of differential equations 
suitable for the Runge-Kutta method. Note that the only 
dimensional quantities in 122|l are the wavenumber k and 
the Hubble function TL. The natural unit for k is there- 
fore iJp. For now, we will ignore the neutrinos, and return 
to them later in section IVII CI 



Initial conditions 



In order to integrate (|22|l numerically, we need some 
initial conditions at the starting time Xi = Ina^, where 
we choose Oi = 10~^. Here we consider only adiabatic 



initial conditions, as derived in pj|: 
<5 = <5h = ?$, 



" 2n 



(23) 



The initial condition for $ acts as a normalization, and 
can be chosen to be 4> = 1 ^. Note that we get 



391 +1^6 = 0. 



(24) 



At early times the optical depth t' is very large, meaning 
that to the lowest order, everything that is multiplied by 
t' in should be zero. This implies 9/ = for I > 2, 
and 9f = for all I. However, when integrating 1)22(1 
numerically, we will need the lowest order non-zero ex- 
pressions for all the multipoles (including polarization). 
And the equations seem to be most well-behaved if we 
also use these very small, but non-zero expressions as 
initial conditions. We therefore derive these expressions 
here (see also [Tslfl^ '). 

Very early, the quantity e = k/iTir') is a small number, 
and can therefore be used as an expansion parameter for 
the multipole hierarchy. As we will see, 9; ~ e9;_i for 



Z > 2, 9^ 



92, 9f - e92, and 9 



e9 



i-i 



for ^ > 3. Assuming that this is true ^, and also that 
the derivatives of the multipoles are of the same order as 
the multipoles themselves, we can compare the order of 
magnitude of the different terms in (|22() . From 9pg and 
9p2 we get the equations 9^ - n/2 = Q2 - n/10 = 0, 
with the result 



9^^ = -9, 



02^=^92, 



H^-92 



Using this in the equation for 92 we then get 

and from the equation for O'pi 

k 



4Ht' 



:92 



(25) 



(26) 



(27) 



^ This does not mean that the perturbation $ of the gravitational 
field has the value 1 (remember that <& is a small quantity), but 
rather that all the other perturbation variables are normalized 
to the value of <I> at a: = Xi. We also ignore the fc-dependence of 
$ at this point, and instead put it back "by hand" in <43l . 

^ This does not lead to any " circular logic" , since we only assume a 
certain asymptotic behavior, and then derive explicit expressions 
proving that the initial assumption was correct. 
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Finally, the equations for Z > 3 reduce to 



we get 



21 + 1 Tir 

I k 



21 



(28) 



(Remember that r' is negative when looking at these ex- 
pressions.) 



D. Tight coupling 

The expressions for the higher temperature and polar- 
ization multipoles in the previous section should be used 
as long as k/{TiT') is small which we refer to as the 
tight coupling regime. However, there's also a more se- 
rious numerical problem in this regime, namely the very 
small value of ZQi+Vb- This quantity is multiplied by t', 
which is very large, meaning that even a tiny numerical 
error in 0i or Vb will result in completely wrong values 
for Q'l and u^, making the system of differential equa- 
tions numerically unstable. This problem can be solved 
by expanding 30 1 + Vb in powers of 1/t', as shown in 
|12lll5l| . Since this is such an important step in order to 
integrate the equations, we include the derivation here. 



Playing around with different parts of 



we get 



r'(3ei + Vb) = 3e; + :^(-eo + 262) - , 

rt n 

(l + RH^-Vb-^^ 

R {3e[+v'b) + ^{-eo + 2Q2)-^'f 



(29) 



3ei + vb = ^ 



ise[ + v',)-v', 

+ ^(-60 + 262)-^* 



(39; + V'f,) + Vb 



(l + i?)r' 



-(-80 + 262) 



(30) 



Taking the derivative of the last equation, using R' = —R 
and substituting various expressions for Vb and w[ so that 
only the combination (30i+Wh) and its derivative appear. 



* We use the tight coupli ng a pproximation as long as Ik/iTir')] < 
1/10 and |r'| > 10 (see llfil 'l. and switch to the full equations no 
later than at the start of recombination (see section fV Bi . 



[(f + i?)r'-i] (3e; + <) 

= - (1 - R:)t'{3Qi + Vb) - (1 + R)t" {3Qi + Vb) 



36'/ 



n 



H'\ k 



n n 



-{-e', + 2e',). 



(-60 + 262) 



(31) 



Until now, everything has been exact. However, during 
tight coupling it is a valid approximation to set 36i + 
equal to zero With x as the variable, this condition 
turns into 



3e'/ + <^-^(36; + «;). 



This gives us the final expression for 3Q'i + u^: 

ise[ + v'b) 



{l + R)r' + ^-l 



(32) 



(33) 



[l - R)t' + [l + R)t" (36i 



Vb) 



n 



1 - ^) ^(-Qo + 2Q2) + ^{-e'o + 26^) . 



This expression is then used in (|29|1 to calculate w^, and 
finally 0[ is obtained from 



1 r 



0'i=o v'b)-v'b 



(34) 



Note that at early times r' ~ 1/a, meaning that r" 



-t'. Therefore, to the leading order, 36^ - 



2(36i 



Vb) 36i +Vb^ a\ 

There is one last technical difficulty in this derivation: 
From H33|) we see that 62 is needed to calculate 36']^ + v'^^ 
and then 6/ . But from (|26|l we see that 6/ is also needed 
to calculate 63. Of course, during tight coupling 62 is 
much smaller than 60, so it is probably a good approx- 
imation to simply set 62 = in l|33|l . Alternatively, one 
could use 62 = as the starting point of a short re- 
currence relation where Q'^ and 62 are calculated with 
growing precision. 



IV. CMB ANISOTROPY SPECTRUM 

When we look at the CMB map today, we are basi- 
cally observing the values of the temperature multipoles 
Qi{k) today. In principle, these can be found by inte- 
grating the system (|22|l of differential equations from Xi 
to X = 0. However, there are two problems that make 
this approach very inefficient. First, we must explicitly 
include all the multipoles up to the highest I we're in- 
terested in, typically I ~ 1200, making the system of 
equations extremely large. Secondly, we must integrate 
the equations for a very large number of values for fc, 
typically several thousand, in order to get an accurate 
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result. This means that even with todays fast comput- 
ers, calculating the CMB anisotropy spectrum would still 
take many hours. Fortunately, the calculation time can 
be reduced by several orders of magnitude by using the 
line-of-sight inte grat ion method, first developed by Seljak 
and Zaldarriaga 



grat i' 



A. Line-of-sight integration 

The basic idea behind the line-of-sight integration 
method is that instead of first expanding H2U|) in mul- 
tipoles and then integrating the equations, we start by 
formally integrating the equation for Q in (|20|) and do 
the multipole expansion at the end. As shown in IJl, we 
first get 



e(fc,/i,77o) = 



')0 



en 



(35) 



Because of the exponential, we can replace ifc/i by d/drj. 
Using partial integrations, expanding H35() in multipoles, 
and using the expression 



M'7-'7o)rf^=^-^[^(^^_,^)] (36) 



for the spherical Bessel functions ji, we then get the fol- 
lowing expression for the multipoles today: 



(37) 



0/(^,??o)= / S{k,T])ji[k{rio-'f])]dri. 
Jo 

The function S{k, rj) is called the source function, 



00 + * + -n 

1 d 



* - $ 



3 



(38) 



ei{k,x = 0)= ' ji[k{r]o - r]{x))]dx 

S{k,x)ji[k{r]Q - f]{x))]dx , (39) 



k drj '^'^ ^ 4fc2 drf ^ 
With X as the variable, this turns into 

r" s{k,x) . 



We therefore need the double derivative of 11 to calcu- 
late the source function. Taking the derivative of the 
appropriate terms in (|22J, we get 



- A. 
dx 

_ 2k 

~ 5n 

3k 
~5H 



QP2 



2k 



5H 
JhC_ 
H 

n 



3fc 

5n 



(42) 
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;e3+ef+er) 



3 






r"n + r'n'" 


+ 10 
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Since we need the derivatives of several perturbation vari- 
ables {Q'i_3,Q'pQ_^,^' and v'^) in order to calculate the 
source function, it can be worthwhile to save these deriva- 
tives along with the variables themselves while integrat- 
ing the differential equations, as these derivatives must 
then be calculated anyway. We are thus avoiding the 
need for methods of numerical derivation. 



B. Calculating Ci 

The observed CMB anisotropy power spectrum today 
is basically given by C/ ~ 'dfix) at the point x — 0, 
i.e. by the Fourier transform of Qf{k). In addition, since 
we have so far ignored the scale-dependence of the initial 
perturbations, we must also include the primordial power 
spectrum P{k). Up to an overall normalization, which 
we ignore for now, the CMB power spectrum is therefore 
given by 

Pik)ef{k). (43) 



Q = 



(27r)3 



With a Harrison-Zel'dovich spectrum predicted by infla- 
tion, the primordial power spectrum is 



27r2 



P{k) 



k 



(44) 



where n is the spectral index, expected to be close (but 
not exactly equal) to 1 from inflation. This gives 



We will return to normalization in section FV El 



(45) 



S{k, x) 



eo + 4' + -n 



1 d 

k dx 



[Ugvb) 



3 d 
4fc2 dx 



V]/' - $' 



The last term in the source function is 



d_ 

dx 



(40) 



&^~gii + mn' (g'n + .gn') 

dx 

+ 71^ (g"n + 2.g'n' + gH") . (41) 



V. PROGRAMMING TECHNIQUES 

We have already mentioned some of the programming 
tricks required to be able to solve all the equations nu- 
merically, including integrating the Peebles equation of 
recombination, and using the tight coupling approxima- 
tion. Here we go through all the other important tech- 
niques needed to get both well-behaved and accurate nu- 
merical solutions, and also to get as fast and efficient 
code as possible. 
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A. Diffusion damping — Boltzmann hierarchy cutoff 

The most obvious thing that has to be done in order 
to integrate is to stop the hierarchy of temperature 
and polarization multipoles at som maximum /max- If we 
choose ^niax large enough, and are careful when selecting 
the cutoff method, there's no need to manually introduce 
a damping scale ~ e~'^' /'^o like the one used in ^J]. The 
easiest cutoff method one can think of is to simply set 
©imax+i — ^Ciax+i ~ Howevcr, this method is poor 
since power is then transferred from /max down to / = 
and back again on a timescale rj ~ /max/fc, because of 
the way the multipoles couple to each other. A very 
high value of /max would therefore be needed to get an 
acceptable result, invalidating the whole purpose of the 
line-of-sight integration method. 

Instead, as discussed in 12], we look at the time de- 
pendence of Qi{k,rj) and Qf{k,rj) for large /, which is 
approximately given by 



(46) 



Now remember the recurrence relation for spherical 
Bessel functions 

21 + 1 

ji+i{x) = ji{x) -ji_i{x) . (47) 

X 

It therefore seems plausible to set 
21 + 1 

e,+i(fc, Tj) ~ ^— e,(fc, 7^) - e,_i(fc, t^) , (48) 

kr] 
21 + 1 

ef+,{k, in) ^ -^^fik, v) - eri(fc, ,7) , / = /max . 

Using this approximation in (|22|l leads to 



l + l 



n 



Hvix) 



(49) 



e'p,^ Aef_^-l±ler + r'er, / = /„ 



n 



nr]{x) 



With this cutoff method, even the low value /max = 6 
gives a good agreement with CMBFAST. Of course, this 
cutoff method is only needed after tight coupling ends, 
since during tight coupling all higher multipoles are ex- 
pressed directly in terms of the lower ones. 



B. Calculating the source function 



much lower resolution after recombination ^. Choosing 
200 points during and 300 points after recombination, 
evenly distributed in x-space, gives a good agreement 
with CMBFAST. It is also sufficient to use 100 different 
values of k between fcmin = O.IHq and fcmax — lOOOi/o 
(for /max — 1200). A bit of trial and error shows that we 
get good results when the /c's are distributed quadrati- 

Cally, that is, ki = fcmin + (fcmax - A:min)(Vl00)2. 

The system of differential equations for each k is in- 
tegrated using an adaptive stepsize fifth-order Runge- 
Kutta method with general Cash-Karp parameters, as 
described in [l3|. We use a relative error of 10~^^. The 
time it takes for the algorithm to finish is roughly pro- 
portional to k, with a maximum of a few seconds for 
k — lOOOi/o- The total time needed to process all 100 
/c- values is therefore about two minutes. This is by far 
the most time-consuming part of the calculation. 

We will later need the source function also at inter- 
mediate values of k and x, that is, we need to make a 
two-dimensional cubic spline. One way of doing this is 
to first take each of the 500 x-values and spline across k. 
Then choose a higher resolution grid of k's, say, 5000 val- 
ues evenly distributed between fcmin and fcmax ^, and for 
each k spline across x. The whole splining process still 
only takes a few seconds to finish. This two-dimensional 
spline is also what requires the most memory in the pro- 
gram - about 120 MB (using 64 bit numbers). 



C. Integrating across x 

The source function is smooth in x, but the Bessel func- 
tion ji in H39|l makes the integrand oscillate for large fc's. 
This may indicate that we should sample the integrand 
at more values of x than the 500 points of the grid. We 
can make a rough estimate of what resolution we should 
use: The Bessel function is really just a combination of 
sin and cos with a period of 2tt. This corresponds to an 
increase in x equal to 

k 

2TT ^ kArj — kri'{x)Ax — ^ Ax , 



Aa; 



2ttH{x) 



(50) 



If we want to have, say, 10 points for each oscillation in 
the Bessel function, we must sample the integrand with 
a resolution 



Aa; = 



2jmix) 
lOfc 



(51) 



The source function S{k, x) in (|40|l is a smooth and 
slowly varying function of both fc and x, except at the 
last scattering surface where it has a sharp peak in x. 
It is therefore sufficient to integrate the system of differ- 
ential equations H22|) and calculate S for a rather small 
number of fc's. For each k the result is stored on a x-grid, 
which has a high resolution during recombination, and a 



^ Here we use the simple definition that recombination "starts" 
when g{x) reaches 10"^" of its maximum value, and "ends" when 
it is reduced to 0.01 of the maximum. For the default model, 
this gives ^start = 1630.4 and z^nd = 614.2. Since g falls of 
exponentially before recombination we don't need to calculate 
the source function before this. 

^ See section fV DI for where the number 5000 comes from. 
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i.e. we must use higher resolution at late times (since 
H{x) is a decreasing function of x) and for large fc's. We 
can also make an estimate of the total number of samples 
in this grid: 



N ■ 



lOfc 



dx 
^^Hix) 



lOfc 
"27 



Vo , Nik 



340iJo) ~ 1800 . 

(52) 



Figure 121 shows a part of the integrand for I — 100 and 
k — 340-ffoj compared to the lower resolution grid with 
500 points in x-space. Clearly the low resolution grid 
fails to sample the oscillations in the Bessel function. 



3 



^ 2 



S 1 
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FIG. 3: The integrand in (I39II plotted as a function of x, for 
I — 100 and k — 340Ho- The top figure shows the entire 
integrand from last scattering until today, using the high res- 
olution X- grid. The bottom figure shows a close-up where the 
oscillations in the Bessel function are important at late times, 
using the high resolution grid (gray line) and low resolution 
grid (black line). In both cases, we have used cubic splining 
between the actual samples of the integrand. 

Amazingly, the CMB power spectrum calculated from 
the low and high resolution grids are indistinguishable. 
This is probably because the dominant part of the x- 
integration comes from recombination, with only a small 
correction from the late time oscillations caused by the 
Bessel function. The source function is non-zero at late 
times mostly due to the integrated Sachs- Wolfe effect, 
which is most important for low Vs (~ low fc's), whereas 
the oscillations in the Bessel function dominate for large 
fc's, hence the details of the oscillations are unimportant. 

The spherical Bessel functions can be calculated us- 
ing algorithms described in [T^ . The argument of the 
function should be chosen from to fcmax'yo ~ 3400, and 



with 10 samples for each oscillation of width 27r we need 
about 5400 samples for each I ^. The result is then splined 
to give a smooth function. The calculation takes a few 
seconds for each I, but since the Bessel functions are in- 
dependent of the cosmological model, they can be calcu- 
lated once and saved to disc for fast access later. 

Finally, the actual integration across x can be done us- 
ing either a simple linear interpolation between the points 
from the (low or high resolution) grid, or using a more 
accurate cubic splining ^. Also in this case, the result- 
ing CMB power spectrum is essentially the same, so we 
therefore choose the faster linear interpolation. 



D. Integrating across k and calculating C; 

The integrand in the final fc-integration, l|45|l . is an os- 
cillating function of k. (This is why we needed the source 
function for more /c's than the ones where we actually in- 
tegrated the differential equations.) Since the dominant 
contribution to the cc-integration is from recombination 
where ?7 ^ 770, we have the rough estimate 

l-rio 

Qi{k) ^ ji{kr]o) S{k,r])dr] const ■ ji{kr]o) , (53) 
Jo 

since the source function varies much slower with k than 
the Bessel function. Thus, 



Ci 



dk. 



(54) 



The Bessel oscillations with period 2tt means that the 
integrand has oscillations with period Afc = 271/770. In 
order to sample each oscillation with 10 points, we must 
therefore use a grid with resolution 



27r 
10% 



(55) 



for the fc-integration. (This leads to the total number of 
fc's (10r/o/27r)(fc,„ax - k^in) ^ 5000 in section lyBl ) 

The computation time can be reduced a bit by ob- 
serving that the full range 0.1 TJq < k < 1000i?o is not 
needed for all I. Instead, from H54|l we see that the peak 
of the integrand is around fc ~ l/rjQ. The integrand falls 
of sharply for smaller fc, but much more smoothly for 
larger fc, as figure 01 shows. As a first estimate, we could 
try the integration range 



0.91 , 21 

< fc < — 

Vo Vo 



(56) 



This estimate 

lOOOiTo and 7?o ~ SAHq-^ (the 
default model). Since we may need larger values of fcmax and 
get larger rjo for other models, it's probably a good idea to use 
at least twice this maximum argument, thus sampling the Bessel 
function at ~ 10 000 points for each I. 

Since we only know the value of the integrand on a grid of finite 
resolution, there's no need to resort to more general methods of 
numerical integration. 
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FIG. 4: The integrand in 11451 plotted as a function of k, for 
; = 100. The peak is slightly after k = //ryo — 29.4J/o, and the 
integrand approaches zero very fast for smaller k. For larger fc, 
however, the integrand dies out much more slowly, and even 
kma.x/Ho ~ I is too Small for an accurate calculation when 
I = 100. With the adaptive algorithm described below, the 
integration continues all the way to fc = 454.9i?o ~ 151/rio. 



Using only this interval gives a rather inaccurate result, 
but we know that this interval at least contains the peak 
of the integrand. One possible algorithm is then to ex- 
tend the interval in steps of one oscillation of width 
Afc = 27r/77o, both to the left and to the right of the 
initial interval, and compare the maximum value of the 
integrand within each step to the global maximum. We 
stop when the local maximum has been reduced to less 
than, say, 10"'' of the global. This algorithm gives a 
CMB power spectrum identical to using the full interval, 
but the fc- integration runs about twice as fast, since on 
average we end up using only half of the interval. 

Since we only do one fc-integration for each I (in con- 
trast to several thousand x-integrations) , we can use the 
slower but more accurate cubic splining instead of linear 
interpolation for the actual integration. The impact on 
the speed of the algorithm from this is negligible. 

Finally, we don't need to calculate C; explicitely for ev- 
ery /. Instead, since the CMB power spectrum is a rather 
smooth function of we only calculate C; for a few Vs 
and use cubic splining to get a smooth function ^. For 
low Vs we should use higher resolution. We choose the 
points I = 2, 3, 4, 6, 8, 10, 12, 15 and 20, and then every 
10th ; up to / = 100, every 20th up to / = 200, every 
25th up to Z = 300, and finally every 50th / above that. 
This gives a total of 44 C;-calculations with Zmax = 1200, 
and the entire integration (across both x and k) only 
takes about 20 seconds with precalculated Bessel func- 
tions. Thus, the calculation of the CMB power spectrum 
is completed in about two and a half minutes. 



E. Normalization 

The final point that has to be considered conserns the 
normalization of the entire power spectrum. The power 
spectrum must be properly normalized if we want to com- 
pare it to CMBFAST. If we want to compare just one 
single model, we could simply use the height of e.g. the 
first peak as normalization. However, since this height 
depends on cosmological parameters, we must be more 
careful if we want to compare several different models 
within the same figure, like in figure [S] IHl and [71 

CMBFAST uses the COBE normahzation which 
basically normalizes to the observed spectrum from 
COBE. The idea is to use a least squares fit of the spec- 
trum for Z < 20 '" (since COBE is only accurate on this 
scale) to a quadratic function in a; = logj^g ^■ 

l{l + l)Ci ~Di[l + D'{x - 1) + D"{x - 1)2/2] . (57) 

From their definition D' and D" are independent of nor- 
malization, and parametrize the shape of the spectrum. 
The fit to the COBE data is then given approximately 
by the formula 

10"Cio = 0.64575 + 0.02282D' + 0.01391(L»')^ (58) 

-o.ol8l9i:'"-o.oo646i:''i:)"-l-o.ool03(i:'")^ 

i.e. the value of Cio is fixed by this expression, and the 
normalization of the rest of the spectrum then follows by 
multiplying the calculated C/'s by the appropriate con- 
stant. One should be aware that this normalization may 
actually introduce a quite significant uncertainty when 
comparing the spectrum to CMBFAST. Part of the rea- 
son is that inaccuracies in the calculation of the low Vs 
get transferred to the entire spectrum by this method. 
By fine-tuning the normalization, the relative error in 
the plots of section IVII can be reduced by up to a fac- 
tor of 2. Since in practice what one really want is to 
fit the calculated spectrum to observational data, and 
not to CMBFAST or some other program, one should 
probably start using the entire spectrum in the normal- 
ization instead of the COBE normalization. Also note 
that CMBFAST gives its output as l{l + l)C//27r, where 
the factor 27r is the commonly used convention. 



^ Note that in the plots of section IVII we do not calculate the 
relative error from the splined C;'s. We use only the explicitly 
calculated C; 's and then cubic splining on the relative error itself. 
This is because the "artificial" error from the splined Cj's is so 
easily removed by simply using more Ts, so this does not really 
indicate any inaccuracy in the algorithms used. 



More precisely, we use the points Z = 3, 4, 6, 8, 12, 15 and 20, 
since these are also the points used by CMBFAST. 
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VI. RESULTS 



Here we compare our program to CMBFAST for some cosmological models. In figureElwe vary the Hubble constant 
h between 0.66 and 0.74, in figurelHlwe vary the baryon density fib between 0.042 and 0.050, in figure[7|we vary the 
CDM density flm between 0.200 and 0.248, and in figure |H1 we vary the spectral index n between 0.95 and 1. In 
figure ini we also plot the power spectrum up to ^ = 2000 for the default model. Following this is the spectrum with 
helium included ffigure lTIil section lVll A|l . a simple model of reionization ffigure lTTl section lVllljll . massless neutrinos 
(figureEl section lVll C|) , and finally the polarization and temperature - polarization cross correlation power spectrum 
(figures US and HI section lyHD)! . 
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FIG. 5: Our program compared to CMBFAST when varying the Hubble constant: the default model h = 0.70 (solid line), 
h = 0.66 (dashed line), and h = 0.74 (dotted line). The figure to the left shows our calculated power spectrum, and the figure 
to the right shows the relative error when compared to CMBFAST (C.F.), which we see is of order 0.3 % or below for most Vs. 
Because the error is so small, there's no point in plotting the spectrum from CMBFAST in the figure to the left, since the 
difference between the curves would be smaller than the width of the lines. 
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FIG. 8: Our program compared to CMBFAST when varying the spectral index: the default model n = 1 (solid line), n = 0.975 
(dashed line), and n = 0.95 (dotted line). 
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FIG. 9: The power spectrum for the default model up to i = 2000, compared to CMBFAST. The relative error starts to increase 
systematically beyond I ~ 1000. In this calculation we have used fcmax = 1500i/o, with 150 values of k chosen initially. We 
have also used /max = 8 instead of 6 in the Boltzmann hierarchy, which reduced the error for I ~ 1800 from 0.9 % down to 
about 0.65%. 
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FIG. 11: The power spectrum with reionization ('section IVIIBfl : Zri — 5 and Az = 0.075 (dashed line), Zri — 10 and — 0.2 
(dotted hne), and no reionization (solid line). We include helium (Yp — 0.24), otherwise the cosmological parameters are as in 
the default model. 
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FIG. 13: The E-mode polarization power spectrum (section IVIIDfl compared to CMBFAST for the default model (solid line), 
n = 0.95 (dashed line), and h — 0.66 (dotted line). The expression for the relative error is slightly modified so it doesn't diverge 
when Ce,i ^ 0: Error = \Ce,i - Cg f | / (C]^ J + e), where 1(1 + l)e/27r = 10"^". The large error of about 3 % close to I = 200 is 
probably due to a slightly wrong position of the minimum. However, the value of Ce,i is still very small here, and the relative 
error doesn't really have a significant meaning until I > 300. 
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FIG. 14; The temperature - polarization cross correlation power spectrum compared to CMBFAST for the default model (solid 
line), n — 0.95 (dashed line), and h = 0.66 (dotted line). Since the relative difference is not very useful for a function that 
crosses zero several times, the figure to the right shows the absolute difference Ac,i = Cc.i — CqJ instead. Note that, with 
the exception of a single point, our program always gives a larger value of Cc,i than CMBFAST. This can not be due to wrong 
normalization, however, since Cc,i is both positive and negative. 



VII. INCLUDING MORE INGREDIENTS 



IS now 



A. Helium 

The primordial mass fraction of ^He, Yp, is defined as 
the ratio between the total mass of helium and the total 
baryon mass. We define a "baryon" as a proton or a 
neutron with a mass mi, ~ mn- Each helium atom thus 
contains 4 baryons, and has a mass approximately equal 
to 4mH, which gives |23| 



M, 



He 



in 



He 



Mb nfcTOH 
nn _ rib - 4nHc 

Ub Hb 



rib 



(59) 



Because of the larger ionization energy, helium recom- 
bines before hydrogen, so that during hydrogen recom- 
bination all helium is essentially neutral. The main ef- 
fect of including helium is therefore that the number of 
free electrons during hydrogen recombination is reduced 
(when keeping Qb fixed). The recombination history of 
helium is well described by the Saha equation. Defining 



Xi 



X2 



"-HC+ + 
"He 



"-H+ 
?1H 



(60) 



we now have three Saha equations [ll 



Xi 



1- Xi - X2 



= 2 



Xi 



2tt J 

3/2 



2tt 



XH _ f meTb\^^^ ^-eo/n 



1 — Xfj V 27r 



(61) 



instead of the one in l|12|l . They are linked in a non- 
trivial way, because the number density of free electrons 



Ue = 2njj(,++ -I- nfjc+ + ?iH+ 

= [2x2 + x,)Yp/ A + xnil - Yp)\ Hb = feUb. (62) 

The ionization energy of neutral and singly ionized he- 
lium is 

Xo = 24.5874 eV , xi = 4 eo = 54.42279 eV . (63) 

Eqs. (|61|) and H62|l are most easily solved by noting that 
they will only be used before hydrogen recombination be- 
comes important, thus fe is of order 1 the whole time 
We therefore use (|61|l to express xi, X2 and a;H in terms 
of fe, and then use H62I) recursively with /g = 1 as a 
starting value. The full machine precision of 15 digits is 
then reached in less than 10 steps. Finally, the electron 
fraction Xg defined in (|10|l is given by 



rie 



fe 



l-Yr, 



(64) 



We switch to the more accurate Peebles equation once 
hydrogen recombination starts fully {X^ < 0.99). At this 
point all helium is neutral, and the only difference from 
section III Bl is that the hydrogen density is now smaller 
than the baryon density. That is, the only changes to 
eqs. (^21 and are 



dXe 

dx H 

nis = (1 - Xe){l - Yp)nb 



f3iTb)il-Xe) - il-Yp)nba^^^iTb)X^, 

(65) 



i.e. nn is replaced by (I — Yp)nb everywhere. The result- 
ing solution Xe is shown in figure IT^ for Yp = 0.24. Note 
that Xe > 1 before helium recombination. 



More precisely, when helium is completely ionized fe = l — Yp/2, 
whereas after helium recombination but before hydrogen recom- 
bination /e ~ 1 — . 
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FIG. 15: The electron fraction Xe. with helium mass fraction 
Yp = 0.24 (solid line) as a function of redshift, compared to 
Yp = Q (dashed line). The first helium recombination He"*"^ — > 
He^ occurs around z ^ 6000, and the second He"^ — > He 
around z ~ 2500. For z > 6000 when helium is doubly ionized 
Xe ~ (1 - Yp/2)/{l - Yp) = 1.1579, and for 2500 <z< 6000 
when helium is singly ionized ~ (1 — 3i^/4)/(l — Yp) — 
1.0789. The residual electron fraction at late times is larger 
by a factor 1/(1 — Yp) with helium included. 



Once the free electron density and the resulting op- 
tical depth have been calculated, the rest of the CMB 
calculation proceeds exactly as without helium. In figure 
El we compare our program to CMBFAST for = 0, 
Yp = 0.24 and Yp = 0.48 (and the other parameters as in 
the default model). The precision with helium included 
is just as good as without. One should also note that the 
other elements (D, ■^Hc, Li etc.) only give corrections to 
the CMB power spectrum of order 10^^ 21] . 



B. Reionization 

At some time long after recombination, we know that 
the hydrogen in the universe became more of less fully 
ionized again. This was probably the result of the en- 
ergetic radiation from the first generation of stars, with 
enough energy to ionize hydrogen, but too low energy to 
ionize helium, which therefore remained neutral. The de- 
tailed mechanism of this process is not fully understood, 
but one possible model is to simply assume that at a cer- 
tain redshift Zri the free electron fraction instantly 
jumps to a constant value, usually = 1, and then 
stays there until today. 

With instant reionization, both the optical depth r' 



and the visibility function g experience a jump disconti- 
nuity at z — Zri, and thus a delta function in g' . This 
can be a bit tricky to implement directly in our program, 
and since it is not very physical either, it is probably 
better to use a smooth (but still sharp) transition from 
the Peebles result to = 1 around Zri- We choose the 
simple formula 



Xe{z) = X, 



Pccblos^^^) . (1 - /) H 

10(z„ - z) 



■ arctan 



Az 



1 

+ 2 



(66) 



where Az can be interpreted as the width of the reioniza- 
tion period, typically chosen to be of order 0.2. Figure 
[TCI shows the free electron fraction, optical depth and 
visibility function with Zri = 10 and Az = 0.2. 
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FIG. 16: The top figure shows the optical depth r (solid 
line) and |r'| (dashed line), and the free electron fraction 
Xe (dotted line), as functions of x. The black lines show 
the semi-instant reionization of eq. 16611 with Zri — 10 and 
Az = 0.2, and the gray lines are without reionization, both 
for the default model with helium {Yp — 0.24). The plateau 
t(z > Zri) — 0.075 is the optical depth to the last scatter- 
ing surface. The bottom figure shows the visibility function 
g — — r'e~^ (solid line), its derivative g'/180 (dashed line), 
and its double derivative g"/65 000 (dotted line). 

Because of the sharp peak in g' , we must use higher 
resolution for the grid of x- values near reionization when 



A more " natural" choice is probably /(^) ~ J e"^'^"^"-'^ dz, hut 
since this function is not entirely trivial to implement in a pro- 
gram, we choose the arctan function instead, which is available 
directly in most programming languages. 
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calculating the source function. We choose to use an ad- 
ditional 200 points between z^^-f Az and z^ — lS.z^ evenly 
distributed in x. The rest of the calculation then pro- 
ceeds as without reionization. The resulting CMB power 
spectrum agrees very well with CMBFAST, even if CMB- 
FAST uses truly instant reionization. We should choose 
Az as small as possible to simulate instant reionization. 
Choosing Az too small, however, leads to problems, since 
the cc-grid must then have an even higher resolution, pos- 
sibly also outside the interval [z^i -I- Az , — Az] . We 
get the best agreement with CMBFAST when Az = 0.2 
for Zri = 10, and Az — 0.075 for z^i = 5. Figure ITTI 
shows the power spectrum compared to CMBFAST. 



C. Massless neutrinos 

The neutrinos decoupled from the cosmic plasma 
slightly before the annihilation of electrons and positrons, 
when the temperature was of order the electron mass. 
The photons were heated by this process, so the neutrino 
temperature is therefore lower by a factor 



1/3 



0.714. 



(67) 



The ratio between the energy densities of neutrinos and 
photons is thus 



P^^^^N -(— 

Pr fir "'SVll 



4/3 



0.681 , (for N^ = 3), (68) 



where N^, is the number of neutrino species, and the fac- 
tor 7/8 is because neutrinos are fermions. Actually, since 
neutrinos are not completely decoupled when the cosmic 
plasma is reheated, one should use an effective number of 
neutrinos Ni, = 3.04 ^2^. With neutrinos as a new com- 
ponent, the Hubble function (0) is of course modified. 
The cosmological constant ft\ in H16|) is also slightly re- 
duced if fif, and flm are fixed. 

The initial conditions for the neutrino monopole and 
dipole are the same as for photons: 



en 



(69) 



The quadropole is more complicated. Since 82 ^ A/2, 



the gravitational potentials are initially related by 

(70) 

0.405 , (for N^ = 3). 



$ = 1^1 
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From 



this gives the initial value 



AA2 = - 



1 



(71) 



Note that A/2 ~ at early times. For the higher mul- 
tipoles we can assume that A/j <C A/;_i. Using l(7T)) and 
H ~ Ho^/^r + f^i/ in 122|l then gives a differential 
equation for A3 that can be integrated directly, with the 
result that A3 ~ . Continuing this way we get A// ~ a' , 
meaning that A/"/ ~ lAfi, and therefore 



{2i + i)n 



l>3. 



(72) 



as the initial condition for the higher multipoles. Finally, 
we use the same cutoff scheme for the neutrino hierarchy 
as for the photons. 



(73) 



We choose ^max = 10 for the neutrino multipoles in order 
to get sufficient precision for large Vs. Note that Hrj ~ 1 
at early times, so that H73I) reduces to A/"/ ~ lAfi directly 
when using (|72|l . 

As a curiosity we find that 



lim A/"'"'* 



(74) 



All the equations where neutrinos appear are therefore 
well-defined in the limit N^, — > 0. The hierarchy of neu- 
trino multipoles is still non- vanishing in this limit, but 
the neutrinos no longer contribute to the CMB power 
spectrum since they decouple from the gravitational po- 
tential. The result is therefore the same as if the neutrino 
hierarchy had not been included at all, as expected. 

The rest of the CMB calculation is exactly the same as 
without neutrinos, as long as one uses the correct expres- 
sion for 5" in the source function (|40|l . Figure IT^ shows 
the power spectrum for = 3, 1 and compared to 
CMBFAST. The agreement is very good below / 600, 
but the error increases somewhat faster for large Vs than 
without neutrinos. 



Here we have implicitly used the fact that photons have two 
polarization degrees of freedom, whereas neutrinos only have one 
(there are no right-handed massless neutrinos). However, the 
neutrino has an antiparticle. The two species thus have the same 
number of degrees of freedom, so there are no additional factors 
of 2 in ISSl . 



D. Polarization power spectrum 

So far, we have only considered the temperature power 
spectrum of the CMB. However, since the radiation is po- 
larized, we also have both a polarization power spectrum 



17 



and a cross correlation power spectrum between temper- 
ature and polarization. With the recent release of the 
three-year results of WMAP , these spectra will be of 
great interest since we now have measurements of the full 
sky CMB polarization map. 

A radiation field in general needs four parameters to 
be described completely, called the Stokes parameters. 
These are the temperature T, linear polarization Q and 
U along two different directions, and circular polariza- 
tion V. T and V are rotationally invariant and can there- 
fore be expanded in spherical harmonics Q and U, on 
the other hand, transform under rotations in the plane 
perpendicular to the direction of the photons. It turns 
out that the linear combination Q ± iU transforms in 
a particularly simple way. Under a rotation tjj it trans- 



where x = k(r]o — rf). Expanding in multipoles, we get 
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forms as (Q ± iU)' 



{Q ± iU), i.e. it has spin ±2 



and can therefore be expanded in what is called spin ±2 
spherical harmonics (see 23] for more details). It also 
means that we can define spin zero quantities by acting 
on Q ±iU twice using the spin raising operator 3 or the 
spin lowering operator 3 (again see [23 for the details) 



Eifi) = -i [g2(Q + iU) + 32(Q - iU) 



\Q + iU)-^''{Q-iU) 



(75) 



The power spectra for E and B are thus rotationally 
invariant, and can be used to describe the polarization of 
the CMB radiation. 

When we are only considering scalar perturbations, we 
can choose a coordinate system (for each Fourier mode) 
where U = Q. We then have Q = Qp and WQ = S^Q 
since 8p only depends on the polar angle. Thus we only 
get E-mode polarization from scalar perturbations. (Ten- 
sor perturbations generate both E- and B-mode polariza- 
tion.) We use the line-of-sight integration method for po- 
larization, similar to the temperature, and get from H2(J|I 



ep(/c,/^,?7o) 
This gives (23] 



1 r^" 

- t(1 - ■P2)ne''=''("-"°)~^dr/ 

2 Jo 

gn(l-/?)e*'=^(''-'"')d77. (76) 



2 

3 



3 rvo 

3 /■"" 

4 70 



JQ 

5n(l + a2)2(x2e-^''")d77, (77) 



Circular polarization can not be generated through Thomson 
scattering, so we will ignore V from now on. 

An alternative method is to construct a 2 X 2 symmetric trace- 
less tensor fro m Q and U, and expand this in tensor spherical 
harmonics |24||. 



ef(fc,%) 



SEik,!]) = 




gU{l + dl)'[x'3i{x)]drj 



gU — —drj 



SE(k,ri)ji[k{rio - r])]d'r] , 



(78) 



Here we have used the result (1 + d^)^[x^ji{x)] = {I — 
l)l{l+l){l+2)ji{x)/x'^, which follows from the differential 
equation j'/ + 2i[/x + [l — l{l + l)/a;2] j/ = satisfied by 
the spherical Bessel function. The E-mode polarization 
power spectrum and its cross correlation with tempera- 
ture is then finally given by [2^ 



C 



E,l 



k 



n-l 



dk 



r°° / k \ """^ dh 

Cc, = I e,{k)ef{k)^. (79) 

Figure [TSl shows the E-mode polarization and figure [T^ 
the temperature - polarization cross correlation com- 
pared to CMBFAST for a few models. The calculation 
uses exactly the same algorithms and techniques as for 
the temperature, only with /max — 8 instead of 6 in the 
Boltzmann hierarchy to get acceptable precision for the 
polarization multipoles. 



VIII. CONCLUSION 

We have here presented all the main steps required in 
writing a program that calculates the CMB anisotropy 
power spectrum. Our focus has been on the computer- 
technical side of the problem, by including all the small 
details that make the program actually work, something 
which is often left out in the literature. We have consen- 
trated on the ACDM model, where the program achieves 
an accuracy comparable to CMBFAST ^'^ (-- 0.1% - 
0.4 %) over a range of cosmological parameters. The pro- 
gram runs in a couple of minutes on a mid-range personal 
computer (as of 2006). While certainly not as good as 
CMBFAST, this is still acceptable considering that the 
code hasn't really been optimized for speed. 

The purpose of this work has been to give a running 
start to those needing to calculate the CMB power spec- 
trum for some exotic cosmological model where the stan- 
dard programs can't be used. With the growing precision 



^'^ See |2j for the details on the extra factor ^(Z - 2)!/(/ - 
I'' The accuracy of CMBFAST is of order 0.1% I2I. 



'2)! 
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of the observed spectrum, a calculation to within the 1 % 
level is often what distinguishes the models and makes it 
possible to rule out some of them. We hope this work will 
encourage others by showing that writing a program from 
scratch to within this accuracy is not really as difficult 
or time-consuming as one may think. 
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